clear all

nreps = 1000;

bootgammatsp  = zeros(2 ,nreps);
bootgammatsp_rob1  = zeros(1 ,nreps);
bootgammatsp_rob2  = zeros(1 ,nreps);


save tspbres bootgammatsp bootgammatsp_rob1 bootgammatsp_rob2

for b = 1:nreps
b;
%set seed
seed = 45454 + 10*b;
rand('state', seed);



% TSP results
load tsp_data


num_coef = 3;

tempindex = find(year >1970 & year < 1987);
cons = ones(length(tempindex),1);
[~,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
temp = randsample(r,length(r),'true');
[temp,~,r,~,~]= regress(tsp(tempindex)+temp, [cons L_tsp(tempindex) year(tempindex) ]);
arvals(1) = temp(2);
temp = sqrt(diag(inv(length(cons))*(r'*r)*inv([cons L_tsp(tempindex) year(tempindex) ]'*[cons L_tsp(tempindex) year(tempindex) ])));
arsevals(1) = temp(2);


tempindex = find(year < 1971);
cons = ones(length(tempindex),1);
[~,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
temp = randsample(r,length(r),'true');
[temp,~,r,~,~]= regress(tsp(tempindex)+temp, [cons L_tsp(tempindex) year(tempindex) ]);
arvals(2) = temp(2);
temp = sqrt(diag(inv(length(cons))*(r'*r)*inv([cons L_tsp(tempindex) year(tempindex) ]'*[cons L_tsp(tempindex) year(tempindex) ])));
arsevals(2) = temp(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = 7; beta = .95; M= 1;

R = length(arvals);
pxest = zeros(2,1);
gammatsp  = zeros(R,1);
dgammatspdrho = zeros(R,1);

for i = 1:R
    pxest(2) = arvals(i);
    [gammatsp(i), dgammatspdrho(i)] = gammafunc_se(T,beta,M,pxest);
end



gammatspse = arsevals'.*dgammatspdrho; 

alphaf = 18.23./gammatsp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Robustness Check 1
tempindex = find(year >1970);
cons = ones(length(tempindex),1);
[~,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
temp = randsample(r,length(r),'true');
[temp,~,r,~,~]= regress(tsp(tempindex)+temp, [cons L_tsp(tempindex) year(tempindex) ]);
arval = temp(2);
pxest = zeros(2,1); pxest(2) = arval;
gammatsp_rob1 = gammafunc_se(T,beta,M,pxest);

% Robustness Check 2
load tsp_data_rob


tempindex = find(year >1970 & year < 1987);
cons = ones(length(tempindex),1);
[~,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) LL_tsp(tempindex) ]);
temp = randsample(r,length(r),'true');
[temp,~,~,~,~]= regress(tsp(tempindex)+temp, [cons L_tsp(tempindex) year(tempindex) LL_tsp(tempindex)]);
ar2val = temp;
pxest = zeros(4,1); 
pxest(2) = ar2val(2);
pxest(4) = ar2val(4);
gammatsp_rob2  = gammafuncAR2(T,beta,M,pxest);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load tspbres

bootgammatsp(:,b) = gammatsp;
bootgammatsp_rob1(b) = gammatsp_rob1; 
bootgammatsp_rob2(b) = gammatsp_rob2; 

save tspbres bootgammatsp bootgammatsp_rob1 bootgammatsp_rob2

end

segammatsp = std(bootgammatsp'); 
segammatsp_rob1 = std(bootgammatsp_rob1);
segammatsp_rob2 = std(bootgammatsp_rob2);

save tspbres bootgammatsp bootgammatsp_rob1 bootgammatsp_rob2 segammatsp segammatsp_rob1 segammatsp_rob2